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was directed to “Numerical Simulation of Shock-Induced Combustion Using Shock-Fitting 
Technique.” Important results of this study were presented at the 30th AIAA/ASME/SAE/ASEE 
Joint Propulsion Conference, Indianapolis, IN, June 27-29, 1994; ALAA Paper No. 94—3100, 
June 1994 (see Appendix A). 
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INTRODUCTION 


One of the propulsion concepts for hypersonic vehicles is the proposed Oblique Detonation 
Wave Engine ([1-3]). In this concept, an oblique shock wave in the combustor is employed 
to increase the temperature of premixed fuel and air to a point where chemical reaction can 
start. It can be utilized to provide a smaller, light weight engine or to provide a higher payload 
capability for a given vehicle weight. 

In the past, many researchers have conducted ballistic range experiments to study supersonic 
combustion/detonation. In these experiments, projectiles were fired in different fuel-air mixtures, 
and detonation structures around the projectiles were recorded. Every gas mixture has a 
detonation wave velocity known as Chapman-Jouget (C-J) velocity, which is characteristic of 
the mixture. If the velocity of the projectile is above the C-J velocity of the reactive mixture, 
the free-stream velocity is referred to as overdriven. On the other hand, if the projectile velocity 
is lower than the C-J velocity, the free-stream velocity is referred to as underdriven. The 
detonation wave structure is highly unstable for projectile velocities less than the C-J velocity 
of the mixture and for a particular projectile diameter. If the projectile is flying above the C-J 
velocity of the gas mixture, the detonation or reaction front structure shows a coupled shock- 
deflagration system near the stagnation line of the body. These two fronts separate from each 
other as one moves away from the stagnation line. The separation between the two fronts occurs 
as soon as the velocity component normal to the bow shock is equal to the detonation velocity. 
The separation between the bow shock and the reaction front is called the induction zone. In 
1972 Lehr [4] conducted a detailed experimental study for a wide range of projectile shapes 
and combustible mixtures. The projectile shapes tested included not only spheres but cones, bi- 
cones, and flat-nose projectiles. The mixtures included hydrogen-air, hydrogen-oxygen, methane- 
air, and methane-oxygen. Ballistic range shadowgraph pictures for Mach 5.11 and Mach 6.46 
from Lehr’s experiments are shown in Figs.l and 2, respectively. In both cases, a free-stream 
temperature of 292 K and a pressure of 42663.2 N/m 2 (320 mm of Hg) is used along with a 
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stoichiometric mixture of hydrogen-air. The projectile diameter was 15 mm. At these conditions 
the C-J Mach number of the mixture is 5.11. Figure 1 shows two discontinuities separated 
from each other. The outer front is the bow shock and the inner front is the reaction front 
produced by ignition of the heated F^-air mixture. The separation between the two is minimum 
near the stagnation point and increases as the shock curves around the body, due to increase 
in induction distance (decrease in post shock temperature) away from the stagnation zone. A 
close examination of the shadowgraphs reveals that as the flow crosses the bow shock, the color 
changes from light to dark, indicating an increase in density. But, as the flow crosses the reaction 
front, the color changes from dark to light, indicating a decrease in density across the reaction 
front. This is due to a large release of energy across the reaction front, causing an increase in the 
temperature; since the pressure remains relatively constant, the density must decrease. Another 
interesting feature is the presence of corrugation in the reaction front. These corrugations are 
caused by the pulsation of the reaction front. The frequency of this pulsation was determined 
to be 1.96 MHz [5]. Figure 2 is for the Mach 6.46 case, and it is seen that the reaction front is 
coupled with the shock near the stagnation line and up to about 60-65 degree body angle from 
the stagnation line. This is because of a very high post-shock temperature at Mach 6.46 that 
causes the induction zone to become so narrow that it appears that the two fronts are merged with 
each other. Decoupling begins further downstream from the stagnation line when the post-shock 
temperature starts decreasing and, therefore, the induction distance increases. Further, both the 
bow shock and the reaction fronts are smooth without any visible instabilities. Thus for an 
overdriven case of Mach 6.46, the instabilities have disappeared. References [4,5] show other 
underdriven cases where it has been shown that the instabilities in the reaction front become 
more pronounced as we reduce the projectile velocity lower than the C-J velocity of the mixture. 
In all these cases the projectile diameter was fixed as 15 mm. 

McVey and Toong [6] also conducted similar experiments where projectiles were fired 
into lean acetylene-oxygen and stoichiometric hydrogen-air mixtures. They proposed the wave 
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interaction model to explain the instabilities in the structure of the detonation wave. Their 
model explains how compression waves can be formed when a new reaction front develops in 
the induction zone between the normal segment of the bow shock and the original reaction front. 
These compression waves lead to a cyclic process which is compatible with most of the observed 
features of the flow. However, the strength of the compression waves remained unresolved in 
their wave-interaction model, which is an important factor in determining if such a model is 
physically possible. Alpert and Toong [7] included the effect of the strength of the compression 
waves and proposed a modified form of the wave-interaction model. 

Several researchers [8-20] have attempted to numerically simulate Lehr’s ballistic range 
experiments [4]. Wilson et al. [8] conducted a detailed numerical investigation of the shock- 
induced combustion phenomena. They used Euler equations and a 13-species and 33-reactions 
chemistry model. They showed the validity of the reaction models and the importance of grid 
resolution needed to properly model the flow physics. Highly resolved calculations for Lehr’s 
Mach 5.1 1 and Mach 6.46 cases with adaptive grids were also performed. The calculations were 
not time accurate and, therefore, the unsteady behavior was not captured. However, for cases 
lower than Mach 5.11, they could successfully capture the instabilities. 

Sussman et al. [9-10] also studied the instabilities in the reaction front for a Mach number 
of 4.79. They also used Euler equations and a 13-species and 33-reactions chemistry model. 
They proposed a new formulation based on logarithmic transformation, which greatly reduces 
the number of grid points needed to properly resolve the reaction front. The unsteady case was 
successfully simulated; however, the frequency was slightly underpredicted. 

Matsuo and Fujiwara [11-12] have studied the instabilities of shock-induced combustion 
around an axisymmetric blunt body. They used Euler equations and a simplified two-step 
chemistry model. They investigated the growth of periodic instabilities by a series of simulations 
with various tip radii and showed that these periodic instabilities are related to shock-standoff 
distance and induction length. They proposed a new model based on McVey and Toong’ s model 
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[6]. The instabilities in the reaction front were explained by their model. 

Ahuja et al. [13-14] and Singh et al. [15] used the Navier-Stokes equations with 7-species 
and 7-reactions [13-14] and 9-species and 18-reactions [15] H 2 -air reaction model to simulate 
Lehr’s Mach 5.11 and 6.46 cases using the shock-capturing method. The Mach 5.11 case was 
found to be unsteady while the Mach 6.46 case was macroscopically stable. The frequency of 
oscillations was found to be in good agreement with the experimental frequency. 

The key parameters for the triggering of instabilities has been identified by various parametric 
studies [16]-[18]. Matsuo and Fujiwara [16] and Ahuja and Tiwari [17] showed that an 
underdriven case, which shows instabilities in the reaction front, can be made stable by having an 
appropriately small size projectile and an overdriven case can be made unstable by having a large 
size projectile. Kumar and Singh [18] concluded that the key parameters for triggering these 
instabilities were projectile velocity, activation temperature, projectile nose radius, reaction rate 
constant, and heat release. 

Tivanov and Rom [19] conducted an analytical study based on an energy equation and a 
chemical rate equation for the flow of a detonable gas mixture over a blunt body. They evaluated 
the conditions for the stability of the detonation process and the appearance of the oscillations. 
The frequency of oscillations matched very well with the experimental data. 

Matsuo et al. [20] simulated the regular and large disturbance regime cases of Reugg and 
Dorsey [21] using two-step chemistry model. With a series of simulations the large disturbance 
regime was explained with Alpert and Toong’s model. Their results revealed that the intensity 
of heat release was a key parameter in determining the regime of the unsteady flow. Flow 
features of the unsteady combustion with low-frequency and high-amplitude oscillation, known 
as large-disturbance regimes, are reproduced when the concentration of the heat release is very 
high. For moderately high heat release, a high-frequency, low-amplitude periodic unsteadiness 
that belongs to regular regimes was observed. 

In problems like shock-induced combustions where physical instabilities are present, the 
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shock-capturing methods will smear some of the instabilities. Thus shock-capturing methods, 
when used in complicated problems of practical interest, would not reproduce many of the 
intricate flow features. On the other hand, in the shock-fitting approach, one knows the precise 
location of the discontinuity which acts as one of the boundaries of the flow field, across which 
Rankine-Hugoniot conditions are applied. This approach ayoids taking differences across the 
shock and the smearing of the shock that occur in the shock-capturing method. There are some 
obvious advantages of shock fitting over shock capturing. Shock fitting requires far less grid 
points compared to shock capturing. In shock-capturing the bow shock becomes a smeared shock 
surface and, requires more grid points for the extension of the grid in the free-stream region. 
This adds to the savings in computational time in shock fitting as compared to shock capturing. 
Since very small dissipation is needed in shock fitting, the intricate details of the flow can be 
reproduced, as the dissipation does not smear the important flow features. 

The present study investigates, in detail, the shock-induced combustion phenomena for 
a premixed stoichiometric H 2 -air mixture flow at hypersonic speeds (Mach 5.11 and 6.46), 
using the shock-fitting technique past a 15 mm spherical projectile. These are the first such 
simulations done with shock fitting for the ballistic range experiments. The analysis is carried out 
using the axisymmetric version of the SPARK2D code [22] with shock-fitting capability, which 
incorporates a seven-species, seven-reactions combustion model for hydrogen-air mixtures. 

BASIC GOVERNING EQUATIONS 

The physical model for analyzing the flow field is described by the Navier-Stokes and species 
continuity equations. For two-dimensional axisymmetric flows, these equations are expressed 
in physical coordinates as 


dU OF dG 

dt dx dy 


( 1 ) 
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where vectors U, F, G, and H are written as 
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The coefficients A;, A, A, A, and A for each species is found by a curve fit of the data 
tabulated in reference 26. 

In Eq. (1) only (N s — 1) species equations need to be considered in the formulation since 

the mass fraction of the species is prescribed by satisfying the constraint equation 

N s 

E/' = 1 (12) 

i = 1 

The specific heat at constant pressure for each species is prescribed in Eq. (11) by a fourth- 
order polynomial in temperature. The binary diffusion equation for the diffusion velocity of 
the i th species 


ui = u{i + ujj 


(13) 
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is as follows: 
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It may be noted that this equation has to be applied only to (N s — 1) species. The diffusion velocity 

N s 

for the remaining species is prescribed by satisfying the constraint equation J2 fwi = which 

?=i 

ensures the consistency. In Eq. (14), it has been assumed that the body force vector per unit 
mass is negligible. In addition, thermal diffusion is considered to be negligible when compared 
with the binary diffusion coefficient. 


CHEMISTRY AND THERMODYNAMIC MODELS 

Chemical reaction rate expressions are usually determined by summing the contributions 
from each relevant reaction path to obtain the total rate of change of each species. Each path is 
governed by a law of mass action expression in which the rate constants can be determined from 
a temperature-dependent Arrhenius expression. In vector H, the term cv t = MiC l represents the 
net rate of production of species i in all chemical reactions and is modeled as follows : 
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where Eq.(15) is a representation of an N r -step chemical reaction, and Eq.(16) is the production 
rate for the i th species as determined from the law of mass action. The reaction constants acq 
and Ac b j are calculated from the following equations 
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(19) 


The equilibrium constant appearing in Eq.(18)is given by 
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The forward rate for each reaction is determined from Eq.(17) which is based on the Arrhenius 
law. The appropriate constants Aj,ctj, and ej for the H 2 — air reaction system can be found in 
[23]. The reverse rate is then calculated from Eq.(18). 

The hydrogen-air combustion mechanism used in this work is based on the Jachimowski 
hydrogen-air model [23] which uses seven species and seven reactions. The species are N 2 , O 2 , 
H 2 , OH, H 2 O, O, and H. Each of the seven reactions can proceed in the forward and backward 
directions. The reactions are 

1) 0 2 + H 2 ^OH + OH 

2) 0 2 + H^ OH + O 

3) H 2 + OH ^ H 2 0 + H 

4) H 2 + O ^ OH +H 

5) OH + OH ^ H 2 0 + O 

6) OH + H + M ^ H 2 0 + M 

7) H + H + M^H 2 +M 

The stoichiometric chemical reaction for a hydrogen-air system can be written as 
2H 2 + 0 2 + 3.76N 2 2H 2 0 + 3.76N 2 

For a blunt body moving through a reactive mixture at hypersonic speeds, the temperature of 


( 21 ) 

( 22 ) 
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the fuel-air mixture after the bow shock is sufficiently high to initiate the reaction. Once the 
ignition starts, chemical energy is released and the reaction front is formed. In the induction 
zone, the temperature and pressure remain relatively constant at the post-shock conditions, while 
the concentrations of radicals build up very rapidly. 

METHOD OF SOLUTION 

The governing equations of motion are transformed from physical space (x,y) to the com- 
putational domain (£,77) by the following relations 


T = t 


f = £(t,x,y) 


(23) 
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The final form of the governing equation in the computational domain with time-dependent 
terms is given by 
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Those terms which add to zero analytically, but numerically are not zero, are referred to as 
GCL (Geometric Conservation Law correction) terms in Eq. 24 above. 

BOUNDARY CONDITIONS 

The flow conditions behind the bow shock are determined by Rankine-Hugoniot relations. 
The shock boundary is allowed to move until it reaches a steady state position. Figure 3 shows 
the coordinate transformation used in the shock-fitting procedure where the following notations 
are employed: Let V\ represent the vector component of the fluid velocity normal to, and 
measured with respect to, the moving shock. Therefore, one may express 

Vi = {[Voo~ U s ] • n]n 

~ ^ ^oo^r * ^ T ^ Ug ' * TL ^ 

Consequently, the magnitude of the shock velocity in the direction normal to the body (i.e., 
along the radius) is given by 
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The derivative which appears in Eq. (26) is evaluated by using the second-order central 


difference formula as 


r n = 
se, 


(r n -r n ^ 
V 5e o+i) / 


(2 < j, < nny - 1) 


(27) 


a U) 2A 9 

At the beginning of the predictor step, the shock wave radial distance is computed from the 
Euler predictor equation 
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Therefore, 

V\ — V\ • n 

( 29 ) 

— f — * 

= v\i r + ui^ 


where vi =component of fluid velocity V\ normal to the body (i.e., along i r direction) and u 1= 
component of the fluid velocity V\ tangential to the body (i.e., along ig direction) Therefore, 



( 30 ) 


and 
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Let u ls = velocity component tangent to the shock (i.e., in s direction) and v ln = velocity 
component normal to the shock (i.e., in n direction). Then by applying shock jump conditions, 
we have 


Pl^ln — P2V2TI 

Since tangential velocity is preserved, then 

u 2s = u ls 


( 32 ) 
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Let V2 = resultant velocity after the shock with respect to shock coordinates. Therefore, 


V2 = v 2n n + u 2s s 

The component of V 2 along i~o (i.e., tangent to the body) is given by 
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Similarly, 
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Pressure behind the shock is obtained from the MacCormack scheme at the predictor level. 
Once the pressure is known behind the shock, the normal component of the flow velocity ahead 
of the shock (measured with respect to shock) can be related to the pressure behind the shock 
by manipulating the oblique shock relations which are 

PlVl = P2V2 
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where Vi and V 2 are resultant velocities. Manipulating these equations gives 
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Equations (28), (36), (37), (39), and (40), when written in the notations of the advanced time 
level in terms of the preceding time level, can be written as 
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Note i is normal to the body and varies from i=l at the surface to i=nnx at the shock. Also 
j is along the body and varies from j=l at the stagnation line to j=nny at the outflow boundary. 

SOLUTION PROCEDURE 

Solution procedures are followed in four steps as described below. 

STEP 1 : Initial Solution 

The initial conditions for this calculation are obtained by using an approximate curve fit for 
the location and shape of the bow shock. Newtonian pressure distribution is used at the body. 
The approximate curve fit of Billig [25] is used to find r s and r s q along the shock. To find 
the initial conditions immediately behind the shock, the radial shock velocity r St is set equal 
to zero and Eqs. (25), (42)-(45) are used. The initial flow conditions on the wall are obtained 
using the known wall temperature in conjunction with the pressure from the Newtonian pressure 
distribution. The initial flow conditions at interior grid points are obtained by assuming a linear 
variation between the flow conditions immediately behind the bow shock and the wall conditions. 

STEP 2 : Predictor Step 

At the beginning of the predictor step, the shock wave radial distance is computed from Eq. 
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(30). The pressure immediately behind the shock yPnnxj) ls com puted using the MacCormack 
scheme 
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Once the pressure behind the bow shock is obtained, - and can be computed from 

the normal shock relations given by Eqs. (42) and (43). Similarly the components of the fluid 
velocity behind the bow shock can be found from the oblique shock relations given by Eqs. 


(44) and (45). The remaining unknown T^ - is calculated using the equation of state. This 
completes the predictor step. 

STEP 3 : Corrector Step 

The corrector step is similar to the predictor step except that the shock wave radial distance 
is evaluated using the modified Euler corrector which is 
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and Eq. (46) is replaced by the MacCormack Corrector scheme in which the usual backward 


difference for dG/dy is replaced by a forward difference given by 
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This completes the corrector step. 
STEP 4 : 


Once the calculation of boundary condition at i=nnx (i.e., after the shock) has been performed 
by the shock fitting method, the predictor or corrector steps are initiated at the interior grid points. 
All other boundary conditions are calculated after the predictor or corrector step is completed 
at all interior grid points. 

The flow conditions along the supersonic outflow boundary (i.e., at j=nny) are determined 
by using a second-order extrapolation of interior grid point data. 
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Along the body surface the boundary conditions of no slip, zero pressure gradient, adiabatic, 
and non-catalytic wall were used. 

RESULTS AND DISCUSSION 

Numerical computations were conducted to simulate Lehr’s [4,5] experimental results. The 
physical and free-stream conditions used in the simulations were: 

Mqo = 5.11 and 6.46 
Poo = 42732N/m 2 
Too = 292K 
r n — 15 mm 

Figure 4 shows the typical grid used in the calculation. For the Mach 5.11 case, calculations 
were carried out on a grid of 101x101. Due to close coupling of bow-shock and reaction 
front (i.e., small induction distance) at high Mach numbers, a finer grid was needed to resolve 
the flow field. Therefore, for Mach 6.46, a grid of 201 x 151 was used with 201 points in the 
circumferential direction and 151 points in the normal directions. 

Figure 5 and 6 show the pressure contour plots for the Mach 5.11 case. The complicated 
wave pattern seen in Fig. 5 can be viewed as made up of two types of compression waves. One 
type of compression wave originates from the reaction front while the other has been reflected 
from the projectile body. The reflected compression wave interacts with the original compression 
wave and, at the point of interaction, two new waves are generated. This reflection produces the 
observed cell structure. The compression wave which moves towards the bow shock, overtakes 
it and causes the bow shock to move forward. Thus the kinks appearing on the bow shock are 
due to some of its structure being distorted by the compression waves. Figure 5 also reveals 
that these pulsations in reaction front are strong near the nose region and dissipate near the 
shoulder regions of the projectile. This fact is further supported by Fig. 7, which shows the 
numerical value of pressure contours along the complete body. From the pressure values given, 
it is shown that the pulsation of compression waves that originate near the stagnation line are 
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the strongest, and as one moves towards the shoulder region of the projectile, the pressure is 
reduced to almost atmospheric pressure. 

Figure 8 shows the pressure distribution along the stagnation line. The pressure increases 
from free-stream pressure to 1.332 e+06 N/m 2 as the flow passes the normal shock along the 
stagnation line. The flow then encounters the pressure wave (see Fig. 6) which further raises the 
pressure to 1.375 e+06 N/m 2 . The pressure then drops to 1.342 e+06 N/m 2 as it passes through 
the expansion wave. This pattern repeats itself as the flow encounters a series of compression 
and expansion waves. 

To help understand the temporal nature of these instabilities, attention is now focused on 
the time history of physical variables along the stagnation line. Figure 9 shows the time history 
plot of water mass-fraction along the stagnation streamline. It shows two discontinuities. The 
outer discontinuity is the bow shock, which shows little kinks in the structure, and the inner 
discontinuity is the reaction front. The highly periodic oscillations in the reaction front that 
originate near the stagnation line and then spread downstream are clearly evident. The bow 
shock is at 0.009224 meter and the projectile surface is at 0.0075 meter (projectile surface is not 
shown in the figure) from the center of the blunt body. As seen from the figure, the frequency 
of oscillation (which is inverse of the time period) can be calculated directly from the plot. This 
frequency is 2.0 MHz, whereas the experimental frequency from Lehr’s ballistic data for Mach 
5.11 is 1.98 MHz. Also, the amplitude of the oscillations of reaction front is 8.0 e-05 meters. 
Figure 10 shows the water mass fraction contours. The reaction front can be seen in this figure. 
The instability is characterized by a regular periodic wave motion having a constant frequency 
and amplitude. The contour plots shown here show the spatial variation of instability at one 
point in time. They also show that instabilities are not restricted to the reaction front only but 
are convected towards the body, thus affecting the entire flow field. The reaction is complete 
near the body where the maximum water mass production can be seen. 

Density contours are shown in Fig. 11. It clearly shows the presence of two discontinuities. 
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The outer discontinuity is the bow shock and the inner discontinuity is the reaction front. These 
two fronts can be seen separated from each other. Moving downstream the induction length 
increases because of lowering of post-shock temperature. The oscillations in the reaction front 
can be seen clearly. It is also seen that the density increases just after the shock and then decreases 
as the flow passes through the reaction front, in agreement with the experimental shadowgraph. 

Figure 12 shows the temperature distribution along the stagnation streamline. Following 
a streamline into the stagnation zone, it is seen that the temperature jumps across the shock 
and then stays constant in the induction zone. Past the induction zone, due to the exothermic 
nature of the H 2 -air reaction mechanism, the temperature increases rapidly reaching almost 1 1 
times (3150 K) the free-stream value. The pulsation in the temperature can be vividly seen 
here. To further compare the experimental data with the numerical data, Fig. 13 shows the 
computed shadowgraph of the Mach 5.11 case. It is seen that the bow shock and the reaction 
front are separated from each other near the stagnation line, and this separation keeps increasing 
downstream of the stagnation line. This is what was observed experimentally also. Also, the bow 
shock is quite smooth with very little waviness but the reaction front clearly shows a periodic 
behavior. The instability is characterized by an almost regular periodic wave motion having a 
constant frequency, similar to that observed experimentally in Lehr’s work. 

By means of time history plots, a comparison of the numerical results with the wave- 
interaction model originally proposed by McVey and Toong [6] and further modified by Matsuo 
and Fujiwara [11], can be made. The essential features of this one-dimensional model are shown 
in Fig. 14. In order to understand how the wave interaction model fits with the numerical results, 
we shall have to consider Figs. (14-18) together. 

Figure 15 shows the time history plot for the pressure along the stagnation streamline and 
Fig. 16 shows the enlarged view showing the actual numerical values of pressure. Figures 
17-18 show the time history plots of density and temperature, respectively, along the stagnation 
streamline. By comparing the actual model shown in Fig. 14 with the x-t diagrams of pressure, 
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density, and temperature shown in Figs. 15-18, it is demonstrated below that the model proposed 
by McVey and Toong fits very well with the present numerical calculations. 

As shown in Fig. 14, a contact discontinuity first approaches the original reaction front. 
The gases are hot on the upstream side of the contact discontinuity and comparatively cold on 
the lower side, as clearly seen in Fig. 18. These hot gases behind the contact discontinuity 
begin to react, generating compression or pressure waves that propagate both upstream and 
downstream, as seen in Fig. 15 and the enlarged view in Fig. 16. The compression wave which 
propagates upstream intersects with the bow shock and produces a contact discontinuity behind 
the bow shock (Figs. 17 and 18). The bow shock is stronger after the interaction and, therefore, 
the gas is hotter on the upstream side of the contact discontinuity. The hot gases behind the 
contact discontinuity reduce the induction time and create a new reaction front, thus generating 
another set of compression waves. At a somewhat later time, the contact discontinuity reaches 
the position of the original reaction front, extinguishing the reaction at this point because no 
more unreacted mixture exists there. The rate of energy release is effectively reduced, which 
generates rarefaction waves as shown in Fig. 16. The reaction front begins to recede because 
of increasing induction time of the colder fluid. The compression wave traveling towards the 
blunt body gets reflected from the body, travels back to the bow shock, and interacts with it at 
about the same time that the most recently created compression wave arrives at the bow shock. 
The compression wave and the reflected compression wave from the body interact with the bow 
shock, thus providing a possible mechanism for the creation of another contact discontinuity, 
i.e., secondary contact discontinuity. The gases, being hotter on the upstream side of the contact 
discontinuity, start burning again generating compression waves; the cycle is then repeated as 
shown in Figs. 15. and 16. Matsuo and Fujiwara also emphasized the importance of considering 
the reflection of the compression wave from the body in their calculations. The compression 
wave reflected from the blunt body may not necessarily be in phase with the compression waves 
created by the new energy release front. Thus, once these reflected waves interact, they cause 
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the flow to be not exactly periodic; however, the pulsating energy release front could still be 
nearly periodic. 

The results for the Mach 6.46 case are now presented. This is a superdetonative case, i.e., the 
projectile velocity is higher than the C-J velocity of the mixture. Figure 19 shows the pressure 
contours as well as the wave pattern similar to Mach 5.11. When compared with Fig. 6, it is 
clear that the frequency of the compression waves moving towards the body and moving towards 
the bow shock is much higher. The bow shock and the reaction front are almost coupled with 
each other. Fig. 20 shows the density contours. For this case, a very small induction distance 
occurs as a result of the post-shock temperature remaining significantly high up to some distance 
downstream of the stagnation zone. Away from the stagnation line, the induction distance is 
increased as a result of the decreasing shock strength and post-shock temperature. The bow shock 
has a very crisp and smooth profile. The reaction front, which is smooth up to about 60-65 
degrees from the nose region, is wrinkled with very small amplitude waves downstream. Figure 
21 shows the pressure distribution along the stagnation streamline. There is a jump in pressure 
after the bow shock and then the pressure drops when the pressure wave meets a rarefaction wave. 
It increases again when it encounters another compression wave. After the energy release front, 
there is another jump in pressure. This pressure wave oscillates between a high value (when it 
encounters a compression wave) to a low value (when it encounter a rarefaction wave). Also, 
when compared with Fig. 8 for Mach 5.11, we see that there are more numerous oscillations 
in pressure for the Mach 6.46 case because of higher frequency compression waves generated. 
Figure 22 shows the temperature distribution along the stagnation streamline. The stagnation 
point temperature is 3550 K. The temperature increases abruptly as the gas encounters the bow 
shock. Immediately after the shock, the temperature stays constant for a short distance and then 
begins to increase due to exothermic reaction. Figure 23 shows the computed shadowgraph for 
density. The bow shock and the reaction front remain coupled with each other up to about 60-65 
degrees from the stagnation streamline, as observed experimentally, and then start decoupling. 
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Also, the reaction front shows slight oscillations of very low amplitude. Figure 24 shows the 
time history plot for water mass fraction. The bow shock is at 0.00884 meters from the center of 
the blunt body. The distance between the bow shock and the reaction front is very small. Also, 
as is clear from the figure, the amplitude of oscillations of the reaction front is 2.5 e-05 meters 
which is quite small as compared with the Mach 5.11 case. The frequency of oscillations can 
be computed directly from this figure and it is found to be 2.85 MHz, which is comparable with 
the earlier work [14]. Thus, Mach 6.46 case is a very high frequency but very low amplitude 
phenomena. Experimental fundamental frequency for the Mach 6.46 case is not available. 

CONCLUSIONS 

A shock-fitting technique has been used to simulate the shock-induced combustion past blunt 
projectiles. In such problems which involve instabilities, the shock-fitting technique gives much 
better resolution of the flow features than the shock-capturing technique. The observed flow 
features have been successfully correlated with the one-dimensional wave-interaction model, 
and the frequency of oscillations has been matched with the experimental data as well as with 
earlier investigations. 
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Figure 1 Shadowgraph of a spherical nose projectile moving at 
Mach 5.11 into a premixed stoichiometric hydrogen-air mixture. 


Figure 2 Shadowgraph of a spherical nose projectile moving at 
Mach 6.46 into a premixed stoichiometric hydrogen-air mixture. 
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Figure 4 Typical grid used in the computation (every 4th point shown) 
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Figure 6 Enlarged pressure contours for Mach 5.11 case 
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Figure 7 Pressure contours showing numerical pressure distribution 
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Figure 8 Pressure distribution for Mach 5.11 along stagnation streamline 
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Figure 9 Time history plot of water mass fraction for Mach 5.11 
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Figure 11 Density contours for Mach 5.11 
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Figure 12 Temperature distribution along stagnation streamline for Mach 5.11 
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Figure 13 Computed shadowgraph for Mach 5.11 
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Figure 14 One dimensional wave interaction model 
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Figure 16 Time history plot of pressure along the stagnation line for 
Mach 5.11 showing numerical values for pressure 
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Figure 21 Pressure distribution along stagnation line for Mach 6.46 
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Figure 22 Temperature distribution along stagnation line for Mach 6.46 
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Figure 24 Time history plot of water mass fraction along stagnation line for Mach 6.46 
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